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A Suydam-unstable circular cylinder of plasma with periodic boundary conditions in the axial 
direction is studied within the approximation of linearized ideal magnetohydrodynamics (MHD). 
The normal mode equations are completely separable, so both the toroidal Fourier harmonic index 
n and the poloidal index m are good quantum numbers. The full spectrum of eigenvalues in the 
range 1 < m < m max is analyzed quantitatively, using asymptotics for large m, numerics for all 
m, and graphics for qualitative understanding. The density of eigenvalues scales like m„ lax as 
m max — > oo. Because finite-m corrections scale as l/m^ ax , their inclusion is essential in order to 
obtain the correct statistics for the distribution of eigenvalues. Near the largest growth rate only a 
single radial eigenmode contributes to the spectrum, so the eigenvalues there depend only on m and 
n as in a two-dimensional system. However, unlike the generic separable two-dimensional system, 
the statistics of the ideal-MHD spectrum departs somewhat from the Poisson distribution, even for 
arbitrarily large m max . This departure from Poissonian statistics may be understood qualitatively 
from the nature of the distribution of rational numbers in the rotational transform profile. 
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I. INTRODUCTION 

The general aim of this paper is to compare and con- 
trast the spectrum of eigenvalues in typical integrable 
wave systems (e.g. waves in a rectangular cavity 0) 
with the spectrum of instabilities in a cylindrical plasma 
within the ideal magnetohydrodynamics (MHD) approx- 
imation. This is a first step in understanding the spec- 
tral problem in the complex three-dimensional geometry 
of the class of magnetic confinement fusion experiments 
known as stellarators 0. 

In ideal MHD the spectrum of the frequencies, u>, of 
normal modes of displacements about a toroidal equilib- 
rium is difficult to characterize mathematically because 
the linearized force operator is not compact [3J. In ad- 
dition to a point (discrete) spectrum of unstable modes 
(u; 2 < 0) there are the Alfven and slow-magnetosonic 
continuous spectra on the stable side of the origin (u> 2 > 
0) and the possibility of dense sets of accumulation points 
on the unstable side. In mathematical spectral theory the 
stable continua and unstable accumulation "continua" Q 
are characterized p| as belonging to the essential spec- 
trum. (For a self-adjoint operator L, the essential spec- 
trum is the set of A-values for which the range of L — X 
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is not a closed set and/or the dimensionality of the null 
space of L — X is infinite.) 

There is experimental evidence that ideal MHD is rel- 
evant in interpreting experimental results 0, QJ , but per- 
haps the greatest virtue of ideal MHD in fusion plasma 
physics is its mathematical tractability as a first-cut 
model for assessing the stability of proposed fusion- 
relevant experiments with complicated geometries in the 
pre-design phase. 

For this purpose a substantial investment in effort has 
been expended on developing numerical matrix eigen- 
value programs, such as the three-dimensional TERPSI- 
CHORE H and CAS3D @ codes. These solve the MHD 
wave equations for perturbations about static equilibria, 
so that the eigenvalue lo 2 is real due to the Hermitic- 
ity (self-adjointness [lb]) of the linearized force and ki- 
netic energy operators. They use finite-element or finite- 
difference methods to convert the infinite-dimensional 
PDE eigenvalue problem to an approximating finite- 
dimensional matrix problem. An alternative approach is 
to use local analysis using the ballooning representation 
and to attempt semiclassical quantization to estimate the 
global spectrum [llj, [l2|, LL3| • 

In order properly to verify the convergence of these 
codes in three-dimensional geometry it is essential to un- 
derstand the nature of the spectrum — if it is quantum- 
chaotic then convergence of individual eigenvalues cannot 
be expected and a statistical description must be used 
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This is perhaps of most importance in understanding 
the spectrum in three-dimensional magnetic confinement 
geometries, in particular the various stellarator exper- 
iments currently running or under construction. These 
devices are called three-dimensional because they possess 
no continuous geometrical symmetries, and thus there is 
no separation of variables to reduce the dimensionality 
of the eigenvalue problem. It has been shown that 
the semiclassical limit (a Hamiltonian ray tracing prob- 
lem) for ballooning instabilities in such geometries may 
be strongly chaotic because there are no ignorable coor- 
dinates in the ray Hamiltonian. 

However, the present paper discusses the opposite 
limit, a system with a sufficient number of symmetries 
to make the ray Hamiltonian integrable and the eigen- 
value problem separable. The geometry is the circular 
cylinder, periodic in the z-direction to make it topolog- 
ically toroidal — we shall refer to the z-direction as the 
toroidal direction and the azimuthal, ^-direction as the 
poloidal direction. The study of this separable system 
will provide a baseline for comparison with the three- 
dimensional toroidal case in future work. The overall 
goal of the paper is to determine if the ideal-MHD spec- 
trum falls within the same universality class as that of 
typical waves in separable geometries or, if not, what 
might cause it to differ. 

Berry and Tabor show that the distribution func- 
tion P(s) for the spacing of adjacent energy levels (suit- 
ably scaled) in a generic separable quantum system with 
more than one degree of freedom is exp(— s), as for a 
Poisson process with levels distributed at random. They 
also show that the spectrum of uncoupled quantum os- 
cillators is nongencric even when the frequency ratios 
are not commensurate, in which case P(s) peaks about 
a nonzero value of s (as also occurs in nonintegrable, 
chaotic systems — the "level repulsion" effect). A more 
surprising departure from the Poisson distribution was 
found by Casati et al. Q for waves in a rectangular box 
with irrational aspect ratio, but the departure was very 
small. Level spacing statistics are discussed also in the 
standard monographs on quantum chaos 0, 0, 0, 0] . 

In contrast with quantum mechanics, where the contin- 
uous spectrum arises from the unboundedness of config- 
uration space, the ideal-MHD essential spectrum arises 
from the unboundedness of Fourier space — there is no 
minimum wavelength. This is an unphysical artifact of 
the ideal MHD model because, in reality, low-frequency 
instabilities with |kj_| much greater than the ion Lar- 
mor radius, a;, cannot exist (where kj^ is the projection 
of the local wavevector into the plane perpendicular to 
the magnetic field B). Indeed, ideal MHD breaks down 
in various ways at large |kj_|, with dissipative and drift 
effects coming into play. 

In this paper we do not attempt to model finitc- 
Larmor-radius stabilization, but instead simply restrict 
the poloidal mode spectrum to m < m max and study 
the scaling of the spectrum at large m max . The nature 



of the dispersion relation is such that the toroidal mode 
numbers n relevant to the spectrum are also restricted. 
In a matrix eigenvalue code such as CAS3D or TERP- 
SICHORE our procedure corresponds to using an arbi- 
trarily fine radial mesh but truncating the toroidal and 
poloidal basis set. 

The eigenvalue equation for a reduced MHD model of 
a stellarator is presented in Sec. [HI We study a plasma 
in which the Suydam criterion |20| for the stability of 
interchange modes is violated, so the number of unstable 
modes is infinite. 

Section IIIII is devoted to developing an understand- 
ing of the dependence (the dispersion relation) of the 
eigenvalues on the radial, poloidal and toroidal mode 
numbers, I, m, and n, respectively. As m and n ap- 
proach infinity, keeping /i = n/m fixed, the growth-rate 
eigenvalues asymptote to a constant, the Suydam growth 
rate, depending only on /i and the radial mode number I. 
We use a combination of perturbation expansion in l/m 
and numerical solution of the eigenvalue equation using 
a new transformation to Schrodinger form that is appli- 
cable over the whole range of m, from 0(1) to oo. This 
gen eralizes the approach of Cheremhykh and Revenchuk 
|2l| . which was limited to the m = oo Suydam eigen- 
value problem. We compare some of the asymptotic re- 
sults in |2lj | with our numerical solutions. Our perturba- 
tion expansion shows that the correction to the Suydam 
limit goes as l/m 2 . Contrary to usual experience [22|| . 
our numerical solutions show that the growth rates do 
not always approach the Suydam values from below as 

771 — * OO. 

In Sec. IIVI we examine the part of the spectrum in- 
volving the most unstable modes, which is essentially 
two-dimensional because only the lowest-order radial 
mode, I = 0, contributes. We relate the considerable 
amount of structure observed in the spectrum to the 
Farey sequences of rational values of the rotational trans- 
form (winding number) of the equilibrium magnetic field. 
Low-order rationals have associated eigenvalue sequences 
giving a regular distribution of eigenvalues locally more 
like the spectrum of a one-dimensional system than a 
two-dimensional one. 

In Sec. [V] we derive the analog of the Weyl formula 
for the average density of states, including an asymptotic 
analysis of the large-Z limit. In Sec. I VII we show level spac- 
ing distributions P{s). Since we are interested in large m 
we first try approximating the eigenvalues by their cor- 
responding asymptotic Suydam limit. This gives a very 
singular distribution with a delta-function-like spike at 
the origin [23| due to the extremely degenerate nature 
of the spectrum in this approximation. By contrast the 
distribution for the exact spectrum has no spike at the 
origin, showing that the small 1 / m 2 corrections break the 
degeneracy sufficiently to completely change the statis- 
tics. 

We examine the statistics for the I = and I = 1 spec- 
tra, both individually and combined (in the low-growth- 
rate region where they overlap) . We have examined suffi- 
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ciently large data sets to show convincingly that the sta- 
tistical distributions are not Poissonian, though that of 
the combined I = and I = 1 spectrum is closest. We also 
split the I — spectrum into two halves to remove over- 
lap of spectra arising from different parts of the plasma. 
These split spectra exhibit a much more dramatic depar- 
ture from Poisson statistics, showing that the ideal-MHD 
interchange spectrum is indeed nongeneric in the sense of 
Berry and Tabor [T^ |. 



II. CHOICE OF MODEL EIGENVALUE 
EQUATION 

The grand context of this paper is the three- 
dimensional linearized ideal MHD problem — to solve, un- 
der appropriate boundary conditions, the equation of mo- 
tion 



(1) 



for small displacements £(r, t) of the MHD fluid about 
a static equilibrium state, where p(r) is the equilibrium 
mass density, r is position, t is time, and F is a Hermi- 
tian linearized force operator ,10| under the inner product 
J d 3 x -F-£ and suitable boundary conditions. (Super- 
script * denotes complex conjugation — we can take £ to 
be complex because all the coefficients in F are real, so 
the real and imaginary parts of £ obey the same equa- 
tion.) 

Most modern magnetic confinement fusion experi- 
ments, in particular tokamaks and stellarators, are 
toroidal. Though not guaranteed for arbitrary three- 
dimensional systems, the equilibrium magnetic field B(r) 
is normally assumed to be integrable in the sense that 
all field lines lie on invariant tori (magnetic surfaces) 
nested about a single closed field line (the magnetic axis). 
Within each toroidal magnetic surface a natural angular 
coordinate system is set up, with the poloidal angle 9 
increasing by 2~k for each circuit around the short way 
and the toroidal angle £ increasing by 2tt for each circuit 
the long way. Each surface is characterized by a mag- 
netic winding number, the rotational transform t, being 
the average poloidal rotation of a field line per toroidal 
circuit, (d6/dQ, over an infinite number of circuits. (In 
tokamak physics the inverse, q = is normally used as 
the rotation number.) 

In this paper we study an effectively circular- 
cylindrical MHD equilibrium, using cylindrical coordi- 
nates such that the magnetic axis coincides with the z- 
axis, made topologically toroidal by periodic boundary 
conditions. Thus z and the toroidal angle £ are related 
through £ = z/Rq, where Rq is the major radius of the 
toroidal plasma being modeled by this cylinder. The 
poloidal angle 9 is the usual geometric cylindrical an- 
gle and the distance r from the magnetic axis labels the 
magnetic surfaces (the equilibrium field being trivially 
integrable in this case). The plasma edge is at r = a. 



In the cylinder there are two ignorable coordinates, 9 
and £, so the components of £ are completely factorizable 
into products of functions of the independent variables 
separately. In particular, we write the r-component as 



r£ r = exp(zm#) exp(— in(,)<p{f) , 



(2) 



where the periodic boundary conditions quantize m and 
n to integers and we choose to work with the stream 
function ip(r) = r£ r (r). 

Since the primary motivation of this paper is stellara- 
tor physics, we use the reduced MHD ordering for large- 
aspect stellarators 0, E); averaging over helical ripple 
to reduce to an equivalent cylindrical problem |25l |26| . 
The universality class should be insensitive to the precise 
choice of model as long as it exhibits the behavior typical 
of MHD instabilities in a cylindrical plasma, specifically 
the existence of interchange instabilities and the occur- 
rence of accumulation points at finite growth rates. 

We nondimensionalize by measuring the radius r in 
units of the minor radius of the plasma column, a, and 
the time t in units of the poloidal Alfven time ta = 
Rq^/JIqP/ Bq, where Bq is the toroidal magnetic field and 
/io is the permeability of free space. Thus lu is in units of 
t^ 1 . Defining A = uJ 2 we seek the spectrum of A- values 
satisfying the scalar equation 



Lip = XMip 



(3) 



under the boundary conditions ip(0) = at the magnetic 
axis and ip(l) = 0, appropriate to a perfectly conducting 
wall at the plasma edge. The operators L and M given 
below are Hermitian under the inner product defined, 
for arbitrary functions / and g satisfying the boundary 
conditions, by 



(f,g)= / r(r)g(r)rdr. 
Jo 



(4) 



The weight factor r in the inner product is a Jacobian 
factor coming from d 3 x — rdrd9dz. 
The operator L is given by 



Id 2 d 

L = —{n — mij r— 

r dr dr 



+ 



rn 



(n — mi) — Ds H (n — mi) 

m 



where the Suydam stability parameter D$ is 



D S = -^p'(r)n'(r) , 



, (5) 



(6) 



with e = a/R <C 1 the inverse aspect ratio, p(r) the 
plasma pressure normalized to unity at r = 0, (3q = 
Z^oPo/Bq the ratio of plasma pressure to magnetic pres- 
sure at the magnetic axis, and fl' the average field line 
curvature. Here 
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where the rotational transform is produced by helical cur- 
rent windings making N 3> 1 turns as £ goes from to 
2tt, fl'(r) giving the averaged field-line curvature. (Note 
that e cancels out in Ds.) We use the notation / = rf'(r) 
for an arbitrary function /, so i = rdt/dr is a measure of 
the magnetic shear and * measures the variation of the 
shear with radius. The term is a measure of the "mag- 
netic hill" that allows pressure energy to be released 
by interchanging field lines, thus driving the interchange 
instability. 

The operator arising from the inertial term in Eq. JTJ> . 



M 



- n2 



1 d d 
r dr dr 



(8) 



is easily seen to be positive definite under the inner prod- 
uct Eq. (@J. 

We observe some differences between Eq. (0 and the 
standard quantum mechanical eigenvalue problem Hip = 
Eip. One is of course the physical interpretation of the 
eigenvalue — in quantum mechanics the eigenvalue E = 
Huj is linear in the frequency because the Schrodinger 
equation is first order in time, whereas our eigenvalue A 
is quadratic in the frequency because it derives from a 
classical equation of motion. 

Another difference is that Eq. is a generalized eigen- 
value equation because M is not the identity operator. 
This is one reason why it is necessary to treat the MHD 
spectrum explicitly rather than simply assume it is in the 
same universality class as standard quantum mechanical 
systems. 

Just as in ordinary eigenvalue problems the eigenvalue 
spectrum for the generalized eigenvalue problem is real, 
and the eigenfunctions ipi have a generalized orthogonal- 
ity property 



(ipi,Mifj) = 5. 



h3 



(9) 



where the normalization has been chosen to make the co- 
efficient of the Kronecker S unity. Here i and j denote 
members of the set {I, m, n}, where I is the radial node 
number and the poloidal and toroidal mode numbers m 
and n, respectively, are defined in Eq. @. The nega- 
tive part of the spectrum, A = — 7 2 < 0, corresponds to 
instabilities growing exponentially with growth rate 7. 

Equation Q is very similar to the normal mode equa- 
tion analyzed in the early work on the interchange growth 
rate in stellar ators by Kulsrud |25| . However, unlike this 
and most other MHD studies we are concerned not with 
finding the highest growth rate, but in characterizing the 
complete set of unstable eigenvalues. 



III. INTERCHANGE SPECTRUM 

In this section we discuss the standard unregularized 
ideal MHD spectrum. It is well known that for A > 
the spectrum consists of the Alfven continuum (the slow- 
magnetosonic continuum being removed in reduced MHD 
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FIG. 1: (a) The rotational transform t(r) = l/q(r) defined 
by Eq. 1111 with to = 0.45, *2 = 0.2 (as for all subsequent 
plots); (b) the magnetic shear parameter i = rt'(r). In (a), 
all distinct rational magnetic surfaces * = /1 = n/m are shown 
for m up to 10. 



|27|). On the unstable side of the spectrum, A < 0, it is 
also known that there is an infinity of eigenvalues pro- 
vided the Suydam interchange instability criterion poj 



t 2 4 



(10) 



is satisfied over some range of r in the interval (0, 1), but 
the details of the spectrum do not appear to have been 
published before. 



A. Profiles 

Interchange instabilities occur only for values of m and 
n such that n — mt vanishes (or at least can be made 
very small |2g) and therefore it is important to know 
something about the function The typical profile 

of *(r) in a stellarator is monotonically increasing in the 
interval [0, a) and we shall assume this to be the case 
here (though it is not always true in modern stellarators). 
For the numerical work in this paper we use a parabolic 
profile 



to + *2r 



(11) 



as illustrated in Fig. QJa). 

Given a rational fraction \x = n^/m^ in the interval 
[*(0),*(a)] (where n^and ?n M are mutually prime) there 
is a unique radius such that t(r^) = fi . Any pair of 
integers (m, n)^^ = (vm^, vn^), v = 1,2,3,... satisfies 
the resonance condition 



(12) 



For example, the set of rationals with 1 < m < 
10 in the interval of t shown in Fig. da) is {/i} = 
{1/2, 5/9, 4/7, 3/5, 5/8}, as shown in the figure. 

To understand the global spectrum we also need to 
know something about the pressure profile. In this paper 
we use a broad pressure profile that is sufficiently flat near 
the magnetic axis that the Suydam instability parameter 
G defined in Eq. 1101) goes to zero at the magnetic axis, 
and for which p' vanishes at the plasma edge 



p(r) = 1 - 6 r 5 + 5 r 6 



(13) 



FIG. 2: (a) The pressure profile p(r), Eq. 1131 . used in this pa- 
per and (b) the Suydam criterion parameter G(r), defined in 
Eq. 1101 (solid line), and the instability threshold 1/4 (dashed 
fine), showing nearly all the plasma is interchange unstable. 



FIG. 3: (a) m — oo eigenfunctions for the I = (solid line), 
I — 1 (short dashes) and I = 2 (short and long dashes) modes 
at fi — 1/2, arbitrary normalization, (b) Growth rates 7 vs. 
resonant * = /1. Dashed lines show approximations Eq. (12 U 
(for I = 0) and Eq. ^ (for I = 1 and 2). 



This profile is shown in Fig. Ufa) and the resulting G- 
profile in Fig.|^b). 



Suydam approximation 



B. High m and n 



In this subsection we choose a particular rational sur- 
face and restrict attention to pairs (m, n) from the 
set {(m, n)^^\v = 1,2,3,...} satisfying the condition 
Eq. 113). 

Denning a scaled radial variable x = m(r — r^) /r^, we 
expand all quantities in inverse powers of m, 



-i L (D 



M 



+ m~ 2 M (2) 



•) (14) 



Also, A = A(°) + : 



2 \( 2 \ and similarly for tp. 



The detailed expressions are given in Appendix [S] . 

We then solve Eq. by equating the LHS to zero 
order by order. At 0(m°), as found by Kulsrud [2^, we 
have the generalized eigenvalue equation 



(15) 



where 



£ (0) 



L (0) 

d 

dx 



A (0) M (0) 



2_2 



A (o))ii + i 2 a; 2 ^A( )-^s(16) 
eta; 



with i and Ds evaluated at r M . For A^ ^ < 0, Eq. |JTSJ| 
can be solved to give a square-integrable eigenfunction 
under the boundary conditions tp^ — > as r — > ±00 
when A^ ) is one of the eigenvalues A M) /, / = 0, 1,2, . . ., 
denoting the number of radial nodes of the eigenfunction 
ip(°) — <pu,l- Note that A M ./ depends only on /i = n/m 
and is otherwise independent of the magnitude of m and 
n. We assume that the y Mi ;(r), when combined with the 
continuum generalized eigenfunctions for A^ ^ > 0, form 
a complete set. 



The leading term in the expansion of the eigenvalue in 
1/ro gives the growth rate in the limit m — > 00, known 
as the Suydam approximation. Restricting attention to 
unstable modes, so that 7 = (—A) 1 / 2 is real, we transform 
Eq. Ijl5(l to the Schrodinger form [2ll | 



d 2 ib 

|^ + Q(r/)^ = , 



where 



?o(??|7>M) = G 



j sech 2 



(17) 



F 2 cosh 2 t] , (18) 



with G = G(r,j) defined as in Eq. 0, T = l/iir,,), 
r\ defined through x = 7sinh^/i(r M ), and ip = 
(coshr?) 1 / 2 (/5(x). [In Ref. OH] Eq. l(T5)l is derived from 
the Fourier transform of Eq. I|15|) . but we can also use 
the real-space version as the equation shares with the 
quantum oscillator the remarkable property of having the 
same general form in both Fourier space and real space.] 
Cheremhykh and Revenchuk |2l| (CR) have made an 
extensive study of the eigenvalues of Eq. Ijl7|l using the 
semiclassical quantization condition 



Q (ri) 1/2 d V = (2Z + 1)tt 



(19) 



which follows from the WKB ansatz tp = 
A(rj) exp ±i J Qq^ 2 drj. CR derive several approxi- 
mations, useful in appropriate limits, improving on the 
earlier result of Kulsrud |25j . In this paper we use two 
of their results to compare with numerical solutions of 
Eq. jnj. The first is E q- ( 4 - 5 ) of EH 



4a 



exp 



2a 



1 

4a 2 



(20) 



which [combining the criteria given in CR's Eqs. (4.4) 
and (4.12)] is applicable when a = (G - 1/4) 1 ' 2 > 1/2. 
The second CR result we use is their Eq. (4.7) 



G- {21 + 1)G 
1 + (4GQ- 



1/2 



(21) 
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applicable when G > T 2 » 1. 

As is seen from Fig. [31 Eq. (|21[) gives a remarkably 
good approximation to the growth rate of the most un- 
stable radial eigenmode, / = 0, and Eq. (|2L)[1 gives a good 
approximation for the higher-Z modes (the semiclassical 
quantization being strictly justifiable only for large 
The growth-rate maxima for each I occur close to the 
maximum of G (and hence T), but not exactly owing to 
the i factor in the definition T = 7/*(r^). 

From Eq. I|21|) we see that, provided the Suydam cri- 
terion G > 1/4 is satisfied, there is an infinity of growth 
rate eigenvalues accumulating exponentially toward the 
origin from above (so the A-values accumulate from be- 
low) in the limit I — ► oo. 

Perhaps less widely appreciated (because m and n are 
normally taken to be fixed) is the fact that there is also a 
point of accumulation of the eigenvalues of Eq. at each 
A^ i as TOmax — ► oo with I fixed. To break the degeneracy 
of A(°) we must proceed further with the expansion in 
1/m. 



2. 1/ 



m corrections 



Proceeding with the expansion Eq. (|14|l . the calcula- 
tion goes through much as in standard time-independent 
quantum perturbation theory [28[ e.g.]. 

The lowest order eigenvalues and eigenfunctions are, as 
found in Sec. lIIIBl A^ ) = A p ,; and ip(°> = ip^i(x), respec- 
tively. The 0(1/ m) correction, \^\ vanishes identically 
from parity considerations — Lp^.i(x) is either an even or 
odd function so its contribution to the matrix elements 
of £W anc l Jif W between ip^ and (p^ is even. On the 
other hand, an d are odd, so A^ = 0. (This 

contrasts with the finite-aspect-ratio toroidal case where 
toroidal coupling of Fourier harmonics of different m to 
form ballooning modes leads to a nonvanishing 1/n cor- 
rection mm.) 

The first nonvanishing correction term is thus 



A' 2 



(aM|£ (2) Im,0 

(m,z|£(%,0</M'|£ (1) M> 



Xfj,,l' — A 



(22) 



where the sum over I' is taken to include an integration 
over the continuum. The operators = LS % ' — A^jMW 
are the higher-order generalizations of C^°>, defined by 
Eq. I |lti|l . The m — oo matrix elements of any operator 
T are defined by 



(fj,, l'\ Flfi, I") = / ip* v {x^ip^v, (x) dx 



(23) 



with the eigenfunctions (p^^x) being normalized so that 
(//, l\M^ \fi, I) = 1. Note that, with the operators L 
and M defined as in Eqs. ((SJ and JSJ , £W is Hermitian 
under the inner product used in Eq. I|23[) only for i = 0. 
However, it can be made Hermitian at arbitrary order by 
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FIG. 4: Growth rates 7 iim 

and 1, found by numerical solution of Eq. fSec. llll Ol : (a) 
n/m = 1/2 (m = 2,4,6,...), and (b) n/m = 11/19 (m = 
19, 38, 57, . . .). At high m the dependence becomes linear, in 
qualitative agreement with Sec. HUB 21 
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FIG. 5: Growth rates vs. radial node number I: (a) found by 
numerical solution of Eq. (Sec. 1111 Gil for (m, n) — (2, 1) 
(diamonds) and (8,4) (triangles), normalized to the n — 1/2, 
infinite-m results (filled boxes); (b) fi = 1/2, infinite-m results 
(points) and asymptotic result Eq. 1201 (dashed line). 



the redefinitions L i— > rL and M i— > rM , which puts the 
eigenvalue equation into Sturm-Liouville form. 

As in quantum mechanics [281 e.g.], if is Hermitian 
the contribution of the second term on the right hand side 
of Eq. I|22|) is always negative for the lowest eigenvalue, 



A 



(0) 



because A,/ — X}"' > 0. However, in ideal MHD 



v(0) 



a positive contribution from the first term usually domi- 
nates and the infinite-m mode is most unstable (23. As 

seen in Fig.Q] this is not always the case: A;f.n is negative 
for [i— 1/2, but positive for /1 = 11/19 = 0.578947 . . .. 

The latter value of /j, is very close to the value giving the 
global maximum Suydam growth rate (see Fig.^J). Thus, 
in the special case studied here and in accordance with 
conventional wisdom, the global maximum interchange 
growth rate occurs at m = 00. Both these results are in- 
tuitively reasonable — the eigenfunctions become increas- 
ingly localized as m — > 00, so the highest growth rate is 
obtained by localizing in the "most unstable" region of 
the plasma, where t ss 11/19. On the other hand, modes 
which localize in "less unstable" regions as m — > 00 can 
achieve a higher growth rate at finite values of m because 
their more extended finite-m eigenfunctions overlap the 
more unstable region and tap into the free energy from 
the pressure gradient in this region. 

Since the eigenvalues approach A M ./ as 1/m 2 , there is 
an infinity of modes in the neighborhood of each X^j in 
the limit TO max —> 00. That is, they are finite-growth-rate 
accumulation points of the complete spectrum. Because 
the rationals fj, are dense on the interval (*(ri), «(ra)), 
where (ri, T2) is the region in which the Suydam insta- 
bility criterion is satisfied, and because X^.i in general 
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depends continuously on /x, the accumulation points A^z 
fill the interval (— 7m ax i0) densely. This is the part of 
the unstable spectrum called the "accumulation contin- 
uum" by Spies and Tataronis 4], though "accumulation 
essential spectrum" might be better terminology mathe- 
matically. 



C. Finite m and n 

In order to calculate arbitrarily high or low-m eigen- 
functions we generalize the transformation in Sec. HUB II 
by the change of variable from r to a new independent 
variable r\ such that 

mt(r) — n = 7sinhrj , (24) 

and a new dependent variable ip(ri) such that 

ip= (icoshr/)" 1 / 2 ^^) , (25) 

so that Eq. © becomes the Schrodinger equation 
Eq. ljT7|. but with Qa replaced by 

Q = G{ri)-~ ^-L - 4? cosh 2 r, 



tanhr? di 1 d 2 i 1 / di , 



2i d?7 2i c??7 2 4i 2 V df] 

Where * = rdt/dr is as defined in previous sections, but 
expressed in terms of r\. 

Differentiating Eq. l(2^|) we find 



dr 7 cosh 77 r 
dr\ i to 



(27) 



Thus, in the large-m limit, equilibrium parameters such 
as G and i are slowly varying functions of 77, e.g. di/ dr) = 
0(l/m) and d 2 i/dn 2 = (9(l/m 2 ). Comparing Eq. 
with Eq. I|18|) we see that, to leading order in 1/to, Q = 
Qo but with i now a slow variable rather than a strict 
constant. 

With the simple form for t, Eq. (|llf). assumed in this 
paper, Eq. l(2^|) is easily inverted to give r(rj), and also 
a cancellation occurs between the terms i'{rj) tanh?y/2i 
and i (r])/2i, so that the exact Q is not much more com- 
plicated than Qo- The eigenvalues in Figs. 0] and were 
computed by integrating Eq. i|17|) with Qq replaced by the 
exact Q and with the appropriate finite boundary con- 
ditions. Low-m results were checked against those from 
an untransformed shooting code. The dashed lines rep- 
resent the results of scans through unquantized, noninte- 
ger values of m to show the smooth, but not necessarily 
monotone, functional dependence of 7 on to 



IV. I = SPECTRUM 

The most unstable modes are those with radial node 
number I = 0. Thus we first consider the set So = 




10 15 20 



0.5 0.55 



FIG. 6: (a) Lattice of quantum numbers on which the part of 
the spectrum between the "ground state" and the threshold 
for the entry of the / = 1 mode is defined, and the unbounded 
contours of constant eigenvalue, (b) The same, in /1 = n/m, 
m space. 
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FIG. 7: Growth-rate eigenvalues of I — modes near the 
maximum growth rate vs. fi = n/m. The ensemble shown is 
for m max = 100. 



{A , m ,„|l < m < TO max ,TO/z min < n < TO^ max }, where 
to and n are integers and /i m i n and ^ max are chosen to 
give the desired range of 7. As we shall be rescaling the 
eigenvalues prior to statistical analysis, it makes no dif- 
ference whether we work with the spectrum of growth 
rates 7 or the eigenvalues A = — 7 2 . However the latter 
choice makes it clearer that the analog of the quantum- 
mechanical ground state is the most rapidly growing 
mode — denoting the maximum growth rate of the I = 
mode by 70, the minimum A is Ao = — 7 2 lax = — 7o- 

The spectrum is defined on the fan-like subset of 
the two-dimensional quantum-number lattice depicted in 
Fig. Eta). Also shown are contours of constant 7 (or 
A), regarded as a continuous function of to and n, which 
are seen more clearly in Fig. Elb) . Here we see a strik- 
ing contrast with more generic systems where the 
constant-eigenvalue contours are segments of topological 
circles enclosing the origin. In the ideal-MHD case the 
contours are topologically hyperbolic, with asymptotes 
radiating from the origin toward infinity. 

An interesting representation of the I = spectrum 
is shown in Fig. A great deal of structure can be dis- 
cerned, determined by the number-theoretic properties of 
the interval of /x depicted. For instance, focusing on the 
low-order rational number 4/7 we define spectral subsets 
S {N/M\4/7) = {A , m ,„|TO = M + 7k,n = N + 4k,k = 
0, 1, 2, ... , [(m max — N)/7]}, where [a;] denotes the largest 
integer < x. 

These spectral sequences all accumulate toward the 
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same Suydam eigenvalue \4/r,o as m max — > oo indepen- 
dently of the choice of M and N. However the rapidity 
of this approach is sensitive to the choice of M/N. For 
instance we see in Fig. the most rapidly converging se- 
quence, 5*0(4/714/7), as a set of points accumulating ver- 
tically from below toward the Suydam eigenvalue. Other 
sequences on either side of 5o(4/7|4/7) approach the 
accumulation point obliquely and much more slowly — 
for m max = 100 they visibly have some distance to go. 
The sequence immediately to the left of 5 (4/7|4/7) is 
S (l/2|4/7), while that to the right is S (3/5|4/7), 1/2 
and 3/5 being the immediate neighbors of 4/7 in the 
Farey sequence 0, p. 300] of order 7 (the first order 
at which 4/7 appears), with the xt-values corresponding 
to So(l/2|4/7) and 5 (3/5|4/7) providing the immediate 
neighbors of 4/7 in each higher-order Farey sequence. 

In discussing the structure of the spectrum it is useful 
to partition 5*o into two subsets, 5*0" and Sq , according 
as the points are to the left or right, respectively, of the 
dashed vertical line shown in Fig. [7| passing through the 
point of maximum growth rate. 

The sequences 5 (l/2|4/7) and 5 (3/5|4/7) accumu- 
late toward A4/70, but slower [0(l/m max )] than does 
5 (4/7|4/7) [0(l/m 2 max ), from Sec. 1111 B 21 . Thus there 
is a gap containing within which 5*o(4/7|4/7) con- 

tributes 0(m max ) points to 5q~, while other sequences 
contribute at most a set of O(l) points. 

Within the gap the spectrum 5q~ is essentially one- 
dimensional, being indexed by the single quantum num- 
ber k. In the full spectrum, So = Sq U Sq , 0(m max ) un- 
related eigenvalues from Sq appear in the gap, making 
the spectrum appear more random and two-dimensional. 



V. WEYL FORMULA 

As discussed in Sec. IIII B 21 the overall maximum 
growth rate for the I = and 1 modes (and, we assume, 
for all /) occurs at m = 00. Thus the threshold value 
when a given mode I first starts contributing to the spec- 
trum is at A = — 7 ; 2 , where 7; is the maximum over ^ of 
I). We denote the corresponding value of /j, by /x;. 

For fixed I and large m max the number of eigenvalues 
Ni(p) in an interval of n/m between /i; and /1 is asymp- 
totically equal to the area in the m, n plane [see Fig.EJa)] 
of the triangle bounded by the lines n = fim, n — fiim 
and m = m max . That is, Ni(ji) ~ - /i;|m^ ax . 

Since contours of constant A (or 7) asymptote to lines 
of constant /1 as m — * 00 we can estimate the number 
of eigenvalues between two values of A (or 7) by invert- 
ing the function A Mi / for /i and substituting this into the 
above expression for Ni(fi). The inverse is double- valued: 
zx = /^(A) > [ii and /u ; ~(A) < Then the number of 
eigenvalues between the ground state A MOi ; and A is ap- 
proximately 

Ff(A)^i| M ± (A)- Mo |mL x - (28) 




FIG. 8: (a) Eigenvalue sequence number iV(7) for the com- 
bined spectral set So U £1 in the case m max = 100. The Weyl 
formula for So, N (7) + No (7)1 is shown dashed, (b) Closeup 
of region containing eigenvalues associated with /j, — 3/5 and 
fi = 5/9. 



The asymptotic dependence of the total spectrum S = 
Sq U S£ U Sf U Sf U . . . is thus 

OO 

iV(A) = ££ivf(A). (29) 
z=o ± 

This is the analog of the Weyl formula 0, p. 258] for 
the integral of the smoothed spectral density ("density 
of states" ) . 

Approximating A M) j = A; + ±(<9 2 A Wi i/<9^)(^ - m) 2 we 

get fif(X) = Mi ± V2(A - k) 1/2 /{& s \ lt „i/0tf) 1/ *- Thus 
there is a square-root singularity at each mode threshold. 

A comparison between the Weyl formula for So and the 
set of points {(7at, N)}, where N is the sequence number 
obtained by sorting the set of I = and / = 1 growth-rate 
eigenvalues from largest to smallest, is shown in Fig.[Sfa), 
showing excellent agreement above the threshold for S%. 
The plotted points may also be regarded as the loca- 
tions of the steps in the "staircase plot" of the piecewise- 
constant integrated density of states function N(j), but 
the scale in this plot is too coarse to resolve the staircase 
structure. 

A finer-scale plot is shown in Fig. Hfb) , in which signif- 
icant deviations from the Weyl curve are seen in the mi- 
crostructure. The range shown in Fig.|SJb) is unusual in 
that it contains two well-defined accumulation sequences 
in close proximity. These are associated with low-order 
values of fi occurring on either side of the growth-rate 
maximum near /1 = 11/19 ~ 0.579 — the sequence as- 
sociated with /x = 5/9 ~ 0.556 is in Sq and the one 
associated with fi = 3/5 = 0.6 is in Sq. There are very 
few eigenvalues associated with high-order rational values 
of fi in the range shown and the two low-order sequences 
present are practically unmixed, either with each other or 
with eigenvalues associated with unrelated higher-order 
rational values of \x. [In fact there is only one such high- 
order mode in the region of the accumulation sequences, 
/i = 51/92, the closest approximant to 5/9 in the set cor- 
responding to 5o(l/2|5/9), which causes the slight jump 
seen in the jj, = 3/5 sequence.] Also, the wide gap con- 
taining no eigenvalues is because the intersection of the 
gaps associated with the two low-order rationals is non- 
empty. 

The spectrum near the marginal stability point, 7 = 0, 
will involve the superposition of many branches of ra- 
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FIG. 9: (a) Nearest-neighbor eigenvalue spacing distribution 
for the approximate spectral set g^uydam ugm g ^ min = 0.5044, 
Mmax = 0.6288, m max = 100 (625 eigenvalues), (b) The same, 
for the corresponding set of accurate eigenvalues Sq. 



dial eigenvalue I . To estimate the asymptotic behavior of 
AT (7) as 7 — * we use the approximate dispersion rela- 
tion Eq. (|20[1 . Taking I to be large we see from Eq. <|20[) 
that the Suydam growth rates 7z(/i) are sharply peaked 
about the location of the maximum, fio, of G(r M ), where 
<j{l-i) = (G — 1/4) 1 / 2 is also a maximum. Thus we can 
expand ct(/x) about fj, 



cr(» 



+ O(( A1 - M0 ) 3 ) , (30) 



where (A//) 2 = — 2er max /er"(/io). To leading order all 
other parameters are evaluated at the maximum point 
fi = 1-J.q. The quadratic correction to a need only be 
retained in the term involving the expansion parameter 
I, so, to leading order, 



7 w 70 exp 



Ttl 



2a(p) 



(31) 



where 70 = 4i(r M0 )(cr max /e) exp(-7r/4<r max - l/4cr^ ax ) 
Solving for /i we find 



M; (7) = Mo ± A/i 



Z 



^max(7) 



1/2 



(32) 



where Z max (7) = (2/7r)er max 111(70/7). Substituting 
Eq. (|32|l in Eq. (|29l) and approximating the sum over I by 
an integral we find the leading order asymptotic behavior 
of the number of eigenvalues to be 



N(rr) 



'max "" max 



111 



To 



(33) 



3ir 7 
which diverges logarithmically as 7 — ► 0. 

VI. NEAREST-NEIGHBOR STATISTICS 

Preparatory to the statistical analysis of eigenvalue 
spacing it is standard practice to rescale, or unfold, the 
eigenvalues so as to make their average separation unity, 
thus making possible the comparison of different systems 
on the same footing. 

We can unfold the spectra by using the Weyl formulae 
above, e.g. for A.; G Sj^ we can define rescaled eigenvalues 



Ef 



<(Ai) 



(34) 



For the set So = Sq U S we can unfold with the com- 
bined Weyl function, ^ ± N a . However, for practical 
purposes we have in this section used empirical least- 
square fits of N(j) to a linear superposition of the basis 
functions (7 max - 7) 1/2 , (7max-7)> (7ma X -7) 3/2 , which 
captures the square-root singularity but avoids having to 
invert 7^/. 

When m m ax is large, the great majority of eigenvalues 
M,m,n are very close to the corresponding m = 00 eigen- 
value with the same [i = n/m, A^,;. Thus one might 
suppose that the statistics of the spectrum are asymp- 
totically the same as those of an ensemble c^ u y dam = 
{A„/ m ,o|l <m< m ma x,m/x min < n < m^ max }. 

In Fig. EJa) we show the distribution of nearest- 
neighbor unfolded eigenvalue spacings for iS uydam , and 
in Fig. |Hfb) that for the set So with the correct finitc- 
m eigenvalues. It is seen that the two distributions are 
radically different — even though low-order rational val- 
ues of fi are rare and the distribution is coarse-grained, 
the high- m approximation induces sufficient extra degen- 
eracy that the Suydam spectrum is dominated by a large, 
but spurious, dclta-function-likc spike at s = 0. (The 
range of /x used in Fig. |5| corresponds to the range of 
I = growth rates above the maximum / = 1 rate, in 
which Sq is the only contributor to the spectrum.) 

The reason why finite-m effects are so important, de- 
spite the smallness of the 0(l/m 2 ) corrections found in 
Sec. IIII B 21 is seen from the Weyl formula, Eq. 
which shows that the average eigenvalue spacing in a 
set containing all values of n/m within the range of in- 
terest scales as TO,7 iax , which is the same order as the 
smallest 0(1/ m 2 ) correction within a set containing only 
n/m — const. Thus in the set of accurate eigenvalues Sq 
there is a strong intermingling of eigenvalues with dif- 
ferent n/m that does not occur in the approximate set 

ciSuydam 

This explains why the nearest-neighbor eigenvalue 
spacing distribution in Fig. Elb) is much closer to the 
Poisson distribution exp(— s) obtained for a random dis- 
tribution of numbers on the real line, and also predicted 
for generic separable systems than that in Fig.^a). 
Nevertheless the set of 625 eigenvalues used in Fig. EJb) 
is too small to say convincingly that the distribution is 
or is not Poissonian, so we need to analyze larger data 
sets to determine how close to generic the ideal-MHD 
spectrum is. 

A cutoff at m max = 1000 gives a set Sq containing 
about 62, 254 eigenvalues in the range between the max- 
imum I — growth rate and the maximum I = 1 growth 
rate. [Note the approximately »n max scaling in the size 
of So, as predicted by the Weyl formula, Eq. (|2"5j) .] In 
Fig. HOf a) we show the nearest-neighbor distribution for 
this set. Close examination of the region near the origin 
reveals no trace of the spike seen in Fig.|^a), not even the 





FIG. 10: (a) Nearest-neighbor eigenvalue spacing distribution 
for the spectral set So using fj, m i D — 0.5044, /i max = 0.6288, 
m max = 1000 (62,254 eigenvalues), (b) The Dyson-Mehta 
spectral rigidity for this set (solid line) compared with that 
for the Poisson process (dashed line). 
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FIG. 11: (a) Nearest-neighbor eigenvalue spacing distribution 
for the spectral set Sq for m max = 1000 (37,932 eigenvalues), 
(b) Nearest-neighbor eigenvalue spacing distribution for the 
spectral set Sq (24,412 eigenvalues). 



tiny spike found by Casati et al. [lj for the spectrum of 
waves in an incommensurate rectangular box. However, 
it is clear that the statistics are not exactly Poissonian. 

In Fig.llU f b ) we show the Dyson-Mehta rigidity param- 
eter A 3 (L) [13, pp. 321-323], defined as the least-squares 
deviation of the unfolded eigenvalue staircase N{E) from 
the best-fitting straight line in an interval of length L. 
Again, the behavior is similar to that for the completely 
random spectrum (Poisson process) in that A3 increases 
linearly with L, but the slope is slightly less than the 
1/15 expected for the Poisson processs. 

In order to understand the departure from Poisson 
statistics better, we show in Fig. ^2 the spacing dis- 
tribution for the corresponding sets Sq and Sq . The 
departure from Poisson statistics is now quite striking. 
This is presumably because of the gaps about low-order 
rational values of fi mentioned in Sec. IIVI which leave 
the 1-dimensional accumulation sequences So(fi\fi) un- 
mixed with other parts of the spectrum, so the spacing 
distribution combines aspects of that for a 1-dimensional 
system (peaked at 1) and that for a generic separable 
2-dimensional system (peaked at 0). 

In Fig. I12f a) we show the spacing distribution for the 
I = 1 spectrum, which is seen to be very much like the I = 

spectrum of Fig. llUf a) in its departure from the Poisson 
distribution. However, we might expect that mixing the 

1 = with the I — 1 spectrum will make the levels appear 
more "random" and Fig. I12f b) confirms that the level 
spacing distribution does indeed become more like the 
exponential expected for a Poisson process. 



FIG. 12: (a) Nearest-neighbor eigenvalue spacing distribution 
for the first 72,500 eigenvalues of the I = 1 spectral set Si, 
m max = 1000. (b) Nearest-neighbor eigenvalue spacing distri- 
bution for the mixed spectral set So U Si over the same range 
of eigenvalues as in (a) (total of 95,000 eigenvalues). 



VII. CONCLUSION 

We have demonstrated that the statistical nature of the 
ideal-MHD interchange spectrum deviates significantly 
from the random Poisson process of generic separable sys- 
tems due to the number-theoretic structure of the eigen- 
value distribution. The similarity between the two level- 
spacing distributions in Fig. 1111 which correspond to two 
different parts of the rotational transform profile, sug- 
gest the possibility that there may nevertheless be some 
universality in the statistics. If so, we have found a new 
universality class. 

The crude regularization used in this paper, simply 
restricting the poloidal mode numbers to m < 77i max , is 
not very physical but corresponds closely to what is done 
in the large three-dimensional eigenvalue codes CAS3D 
and TERPSICHORE 0. Thus, apart from funda- 
mental mathematical interest, the primary motivation 
of this paper has been the numerical analysis of the 
three-dimensional ideal-MHD spectrum as produced by 
these codes. Preliminary results 23] on an interchange- 
unstable stellarator test case show spectra with eigen- 
value separation statistics similar to those of strongly 
quantum chaotic systems. However, the results of the 
present paper indicate that some caution should be taken 
in interpreting ideal-MHD spectra in terms of conven- 
tional quantum chaos theory because of the radically dif- 
ferent nature of the dispersion relation. 

In subsequent work it will be important to examine the 
effect of finite Larmor radius on the spectrum. However, 
this typically makes the problem non-Hermitian and less 
easy to compare with standard quantum chaos theory. 
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APPENDIX A: 1/m 2 CORRECTIONS 



and JSJ 



The coefficients of the expansion Eq. (|14l) are found 
by Taylor expansion of the geometric and equilibrium 
quantities in Eqs. (JSJ 



L (0) 
L (2) 



d . 2 2 ^ 



:2 2 



dx 
d 



dx 
d 



x—i x — —itx — 

ax ax ax ax 
„3 
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dx 
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2i 
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7i* 



dx 
i "t 
~3~ 



dx 



6D S + 5£ s - D s . 
♦ (AD 



M< 2 ) 



dx 2 
da; 2 



+ 1 



dx dx 



2x 



-x — j+x— x- — h 3x 
dx A dx dx 



(A2) 



where, as in the main text, dots denote nondimcnsional 
derivatives, i = r^dijdr etc., and all equilibrium quanti- 
ties are evaluated at r = r„ . 
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